maxpop = max(c(out_model_punkty_adresowe$TOT, out_model_punkty_adresowe$EST_POP))
p1 <- ggplot(data = out_model_punkty_adresowe) +
geom_sf(aes(fill = TOT)) +
scale_fill_gradient2(name = "Liczba ludności", low = "white", high = "red", limits = c(0,maxpop)) +
labs(title = "Ryc.1 Mapa oryginalnych wartości liczby ludności") +
theme_bw() + theme(
axis.text = element_blank(),
axis.ticks = element_blank(),
panel.grid = element_blank())
p2 <- ggplot(data = out_model_punkty_adresowe) +
geom_sf(aes(fill = EST_POP)) +
scale_fill_gradient2(name = "Liczba ludności", low = "white", high = "red", limits = c(0,maxpop)) +
labs(title = "Ryc.2 Mapa estymowanych wartości liczby ludności") +
theme_bw() + theme(
axis.text = element_blank(),
axis.ticks = element_blank(),
panel.grid = element_blank())
print(p1)
print(p2)
p3 <- ggplot(data = out_model_punkty_adresowe) +
geom_sf(aes(fill = RES)) +
scale_fill_gradient2(name = "Błąd estymacji", low = "blue", high = "red", limits = c(min(out_model_punkty_adresowe$RES), max(out_model_punkty_adresowe$RES))) +
labs(title = "Ryc.3 Mapa różnic między wartością obserwowaną liczby ludności a jej estymacją") +
theme_bw() + theme(
axis.text = element_blank(),
axis.ticks = element_blank(),
panel.grid = element_blank())
print(p3)
out_model_punkty_adresowe$RES_reclass = case_when(
out_model_punkty_adresowe$RES < 0 ~ "bledy ujemne",
out_model_punkty_adresowe$RES == 0 ~ "brak bledu",
out_model_punkty_adresowe$RES > 0 ~ "bledy dodatnie",
)
p4 <- ggplot(data = out_model_punkty_adresowe) +
geom_sf(aes(fill = RES_reclass)) +
scale_fill_manual(name = "Typ błędu",
values = c("bledy dodatnie" = "red", "brak bledu" = "white", "bledy ujemne" = "blue")) +
labs(title = "Ryc.4 Mapa zreklasyfikowanych błędów modelu") +
theme_bw() + theme(
axis.text = element_blank(),
axis.ticks = element_blank(),
panel.grid = element_blank())
print(p4)
meanerror = mean(out_model_punkty_adresowe$RES)
rmse = sqrt(mean((out_model_punkty_adresowe$RES)^2))
print(paste("średni błąd estymacji modelu wynosi", round(meanerror,2), "natomiast pierwiastek średniego błędu kwadratowego wynosi", round(rmse,2)))
## [1] "średni błąd estymacji modelu wynosi 3.67 natomiast pierwiastek średniego błędu kwadratowego wynosi 30.44"
Maciejo,
YOU KNOW WHAT TO DO
# wykres rozrzutu
ggplot(pop1km, aes(x=TOT, y=N_BUDYNKI)) + geom_point() +
coord_fixed(ratio = 1) +
xlab("Liczba ludności") +
ylab("Liczba mieszkań") +
geom_smooth() +
geom_abline(color = 'grey') +
theme_bw()
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
model_liniowy <- lm(TOT ~ N_BUDYNKI, data = pop1km)
par(mfrow = c(2, 2))
plot(model_liniowy, which = 1)
plot(model_liniowy, which = 2)
plot(model_liniowy, which = 3)
plot(model_liniowy, which = 4)
mtext("Ryc. 5 Wykresy diagnostyczne modelu liniowego", outer = TRUE, line = -1.5, cex = 1)
W jakim stopniu liczba ludności (TOT) jest wyjaśniana przez liczbę mieszkań?
## Statystyki
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -98.6520 -2.5883 -0.7094 0.0000 -0.5883 174.7735
## Średnia: 3.738246e-16
Wartości odstające można zidentyfikować wykorzystując wykresy diagnostyczne, a w szczególności wykres Residuals vs. Leverage.
p <- autoplot(model_liniowy) + theme_classic()
gridExtra::grid.arrange(grobs = p@plots, top = "Ryc. 7 Wykresy diagnostyczne modelu liniowego")
#dopasowanie modelu
smr_model_lm <-summary(model_liniowy)
c( R2 = smr_model_lm$r.squared, R2_adj = smr_model_lm$adj.r.squared)
## R2 R2_adj
## 0.7019419 0.6990481
~70% zmienności zmiennej TOT jest wyjaśniona przez zmienność zmiennej
N_BUDYNKI.
# Określenie dowolnego przedziału ufności dla współczynników modelu
confint(model_liniowy, level = 0.99)
## 0.5 % 99.5 %
## (Intercept) -7.022728 8.441527
## N_BUDYNKI 3.275627 4.603257
outlierTest(model_liniowy)
## rstudent unadjusted p-value Bonferroni p
## 21 9.742598 3.0248e-16 3.1760e-14
## 25 6.813338 6.7995e-10 7.1395e-08
## 34 -5.176442 1.1389e-06 1.1958e-04
ggplot(pop1km, aes(x=TOT, y=N_BUDYNKI)) +
geom_point() +
geom_point(data = pop1km[c(21, 25, 34),], aes(x=TOT, y=N_BUDYNKI), color = "red", size = 4) +
labs(title = "Ryc.8 Wykres rozrzutu z oznaczonymi punktami o wartościach odstających") +
coord_fixed(ratio = 1) +
xlab("Liczba ludności") +
ylab("Liczba mieszkań") +
geom_smooth() +
geom_abline(color = 'grey') +
theme_bw()
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
meanerror = mean(model_liniowy$residuals)
rmse = sqrt(mean((model_liniowy$residuals)^2))
print(paste("średni błąd estymacji modelu wynosi", meanerror, "natomiast pierwiastek średniego błędu kwadratowego wynosi", round(rmse,2)))
## [1] "średni błąd estymacji modelu wynosi 3.73824648485021e-16 natomiast pierwiastek średniego błędu kwadratowego wynosi 26.25"
PYTANIA Czego dowiadujemy się z wykresów diagnostycznych? Czy model jest dobrze dopasowany? Czy spełnione są założenia regresji liniowej?
# Usunięcie wierszy 21, 25 i 34
pop1km_n <- pop1km[-c(21, 25, 34), ]
# wykres rozrzutu
ggplot(pop1km_n, aes(x=TOT, y=N_BUDYNKI)) + geom_point() +
labs(title = "Ryc.9 Wykres rozrzutu po usunięciu punktów z wartościami odstającymi") +
coord_fixed(ratio = 1) +
xlab("Liczba ludności") +
ylab("Liczba mieszkań") +
geom_smooth() +
geom_abline(color = 'grey') +
theme_bw()
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
# model liniowy
model_liniowy_n <- lm(TOT ~ N_BUDYNKI, data = pop1km_n)
p1 <- autoplot(model_liniowy_n) + theme_classic()
gridExtra::grid.arrange(grobs = p1@plots, top = str_wrap("Ryc. 10 Wykresy diagnostyczne modelu liniowego po usunięciu punktów z wartościami odstającymi"))
#dopasowanie modelu
smr_model_lm_n <-summary(model_liniowy_n)
c( R2 = smr_model_lm_n$r.squared, R2_adj = smr_model_lm_n$adj.r.squared)
## R2 R2_adj
## 0.8914342 0.8903486
~89% zmienności zmiennej TOT jest wyjaśniona przez zmienność zmiennej
N_BUDYNKI.
# Określenie dowolnego przedziału ufności dla współczynników modelu
confint(model_liniowy_n, level = 0.99)
## 0.5 % 99.5 %
## (Intercept) -2.413594 3.198192
## N_BUDYNKI 3.198037 3.843294
outlierTest(model_liniowy_n)
## rstudent unadjusted p-value Bonferroni p
## 65 -6.324628 7.3790e-09 7.5266e-07
## 33 4.742867 7.0854e-06 7.2271e-04
ggplot(pop1km_n, aes(x=TOT, y=N_BUDYNKI)) +
geom_point() +
geom_point(data = pop1km_n[c(65, 33),], aes(x=TOT, y=N_BUDYNKI), color = "red", size = 4) +
labs(title = str_wrap("Ryc.11 Wykres rozrzutu z oznaczonymi punktami o wartościach odstających (po usunięciu punktów)")) +
xlab("Liczba ludności") +
ylab("Liczba mieszkań") +
geom_smooth() +
geom_abline(color = 'grey') +
theme_bw()
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
#WIZUALIZACJA
par(mfrow = c(2, 4))
plot(model_liniowy)
plot(model_liniowy_n)
mtext(str_wrap("Ryc. 12 Wykresy diagnostyczne modelu liniowego (pierwszy wiersz - przed usunięciem, drugi wiersz - po)"), adj = 1, line = 20.5, font = 2, cex = 0.8)
meanerror = mean(model_liniowy_n$residuals)
rmse = sqrt(mean((model_liniowy_n$residuals)^2))
print(paste("średni błąd estymacji modelu wynosi", meanerror, "natomiast pierwiastek średniego błędu kwadratowego wynosi", round(rmse,2)))
## [1] "średni błąd estymacji modelu wynosi 2.30805383656037e-16 natomiast pierwiastek średniego błędu kwadratowego wynosi 9.14"
maxpop1 = max(c(pop1km$TOT, pop1km$EST_POP))
ggplot() +
geom_sf(data = pop1km, aes(fill = TOT)) +
scale_fill_gradient2(name = "Liczba ludności", low = "white", high = "red", limits = c(0,maxpop1)) +
labs(fill = "Liczba ludności", title = "Ryc.13 Mapa oryginalnych wartości liczby ludności") +
theme_bw() + theme(
axis.text = element_blank(),
axis.ticks = element_blank(),
panel.grid = element_blank())
# Przewidywanie EST_POP na podstawie modelu liniowego
pop1km$EST_POP <- predict(model_liniowy, newdata = pop1km)
# Tworzenie mapy z estymowaną liczbą ludności (EST_POP)
ggplot() +
geom_sf(data = pop1km, aes(fill = EST_POP)) +
scale_fill_gradient2(name = "Liczba ludności", low = "white", high = "red", limits = c(0,maxpop1)) +
labs(fill = "Estymowana liczba ludności", title = "Ryc.14 Mapa estymowanych wartości liczby ludności") +
theme_bw() + theme(
axis.text = element_blank(),
axis.ticks = element_blank(),
panel.grid = element_blank())
# Obliczenie reszt
pop1km$Reszty <- pop1km$TOT - pop1km$EST_POP
# Tworzenie mapy reszt
ggplot() +
geom_sf(data = pop1km, aes(fill = Reszty)) +
scale_fill_gradient2(name = "Błąd estymacji", low = "blue", high = "red", limits = c(min(pop1km$Reszty), max(pop1km$Reszty))) +
labs(fill = "Reszty", title = "Ryc.15 Mapa różnic między wartością obserwowaną liczby ludności a jej estymacją") +
theme_bw() + theme(
axis.text = element_blank(),
axis.ticks = element_blank(),
panel.grid = element_blank())
# Przeklasyfikowanie reszt do trzech klas
pop1km$Klasy_reszt = case_when(
pop1km$Reszty < 0 ~ "Niedoszacowane",
pop1km$Reszty == 0 ~ "Brak błędu",
pop1km$Reszty > 0 ~ "Przeszacowane",
)
# Tworzenie mapy przeklasyfikowanych reszt
ggplot() +
geom_sf(data = pop1km, aes(fill = Klasy_reszt)) +
scale_fill_manual(name = "Typ błędu", values = c("Przeszacowane" = "red", "Brak błędu" = "white", "Niedoszacowane" = "blue")) +
labs(title = "Ryc.16 Mapa zreklasyfikowanych błędów modelu") +
theme_bw() + theme(
axis.text = element_blank(),
axis.ticks = element_blank(),
panel.grid = element_blank())
## Statystyki reszt:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -98.6520 -2.5883 -0.7094 0.0000 -0.5883 174.7735
## Średnia reszt: -5.399296e-16
## Odchylenie standardowe reszt: 26.37129
# Histogram reszt
ggplot(pop1km, aes(x=Reszty)) +
geom_histogram(bins=30, fill="blue", color="black") +
labs(title = "Ryc.17 Histogram reszt") +
xlab("Reszty") +
ylab("Częstość") +
theme_minimal()
leaflet(data = pop_100m1, width = "100%", height = "600px") %>%
addTiles() %>%
addPolygons(
data = pop_100m1,
fillColor = ~pal(TOT),
group = "Liczba ludności",
color = "#BDBDC3",
weight = 1,
opacity = 1,
fillOpacity = 1,
highlightOptions = highlightOptions(
weight = 2,
color = "#666",
fillOpacity = 1,
bringToFront = TRUE
),
label = ~paste("Liczba ludności:", TOT),
labelOptions = labelOptions(
style = list("font-weight" = "normal", padding = "3px 8px"),
textsize = "15px",
direction = "auto"
)
) %>%
addLegend(
pal = pal,
values = ~TOT,
opacity = 1,
title = "Liczba ludności",
position = "bottomright"
) %>%
addLayersControl(
baseGroups = c("Podkład mapy"),
overlayGroups = c("Liczba ludności"),
options = layersControlOptions(collapsed = FALSE)
) %>%
setView(lng = mean(st_coordinates(pop_100m1)[,1]),
lat = mean(st_coordinates(pop_100m1)[,2]),
zoom = 12)